tsi_1984-2021_DE

1 TSI I/C/O

1.1 Import TSI images (single location)

The processing pipeline is build upon have a distinct raster for each Time series interpolation case, listing this into a single object and then using lapply to perform functions across each TSI case.

#  base case
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/base_case/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSS.vrt"
tsi_base_raw <- rast(link)

#  case day_16_sigma_8_16_32
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.vrt"
tsi_day_16_sigma_8_16_32_raw <- rast(link)


# case day_16_sigma_8_8_8_16_16_16_64
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_8_8_16_16_16_64/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.vrt"
tsi_day_16_sigma_8_8_8_16_16_16_64_raw <- rast(link)

# case day_16_sigma_8_8_8_16_16_16_32_64
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_8_8_16_16_16_32_64/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.vrt"
tsi_day_16_sigma_8_8_8_16_16_16_32_64_raw <- rast(link)


# All cases are combined into a single list to facilitate using lapply based processing
tsi_list <- list(tsi_base_raw, tsi_day_16_sigma_8_16_32_raw, tsi_day_16_sigma_8_8_8_16_16_16_64_raw, tsi_day_16_sigma_8_8_8_16_16_16_32_64_raw)

cases <- c("base_case", "day_16_sigma_8_16_32", "tsi_day_16_sigma_8_8_8_16_16_16_64", "tsi_day_16_sigma_8_8_8_16_16_16_32_64")
names(tsi_list) <- cases 

tsi_obs_length <- list()
lapply(seq_along(tsi_list),
       function(i) 
       tsi_obs_length[i] <<-  length(names(tsi_list[[i]]))
       )

1.1.1 Import DE shp

link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/DE_shp/vg2500_bld_EPSG3035.shp"

DE <- read_sf(link)

1.1.2 Import CSO data

2 test vis of TSI

create stack for direct location comparison unique index for each TSI case as number of layers differ between them

base_case_obs_count <- readRDS(file = "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/base_case_obs_count.rds")

tsi_stack <- c(tsi_list[[1]]$`20000514_LND05`,
               tsi_list[[2]]$`20000509`,
               tsi_list[[3]]$`20000509`,
               tsi_list[[4]]$`20000509`)

plot(tsi_list[[1]]$`20160807_LND08`)
 plot(vect(DE), add=TRUE)

2.1 static maps

#de <- st_read("/data/Dagobah/dc/deu/vector/DEU.gpkg")

#plot(tsi_list[[2]][[800:808]])

for (i in seq_along(tsi_list)) {
  
  plot(tsi_stack[[i]],
       legend=T, col = viridis(n=100,option="D"),
       range=c(0,10000), 
       ext = ext(DE),
       plg=list(title= paste(names(tsi_list[i]), "\n", names(tsi_stack[[i]]))
       ))
  
   plot(vect(DE), add=TRUE)
       
}

2.2 Interactive map

can be enables if desired, but is computationally quite demanding

val = as.numeric(c(0:10000))
pal = colorNumeric(c("viridis"), val,
                na.color = "grey")


#pal <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(tsi_list[[4]][[675]]), na.color = "grey")

leaflet() |>
  addTiles() |>
  addRasterImage(tsi_list[[4]][[675]], group = "show_layer",
                 maxBytes = (20 * 1024 * 1024)) |> 
  addLegend(pal = pal, 
            values = values(tsi_list[[4]][[600:600]]),
            title = "Reflectance") |> 
  addLayersControl(
    overlayGroups = c("show_layer"),
    options = layersControlOptions(collapsed = FALSE)
    )

## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA

3 data wrangeling of TSI images

3.1 Compute and set min and max values for tsi cases

The empty layers are included in the force output to satisfy the specified day window intervals for the interpolation. Unfortunately this has the downstream effect that when attempting to convert the spatial object to a df, R throws errors when empty layers are included. The current workaround is to compute the the min/max per layer, create an index of valid layers, and then use this index to subset the tsi stacks. It is very slow, but slightly faster than any other method tried so far.

lapply(seq_along(tsi_list),
       function(i) 
      tsi_list[[i]] <<-  setMinMax(tsi_list[[i]])
)

# 
# lapply(1:length(tsi_list), function(i)
# 
  # terra::writeRaster(tsi_list[[i]],
  #                    filename= paste0("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/postproccesed_tifs/",
  #                                     names(tsi_list[i]),"_postprocssed.tif"),
  #                    filetype='GTiff',
  #                    overwrite=TRUE)
# 
# )

NA layers currently are only present in the TSI case with the smallest Sigma values (8_16_32).

3.2 remove all NA layers

# create filter out all layers that dont contain any values
tsi_na_index <- list()
lapply(seq_along(tsi_list),
       function(i) 
         
         tsi_na_index[i] <<- minmax(tsi_list[[i]]) |> 
         as.data.frame() |> 
         slice(1) |> 
         pivot_longer(everything()) |> 
         mutate(id = row_number()) |> 
         na.omit()
)

# convert into index vector
lapply(seq_along(tsi_list),
       function(i) 
         tsi_na_index[[i]] <<- as.vector(tsi_na_index[[i]])
       
)

#saveRDS(tsi_na_index, file = "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/tsi_na_index.rds")
tsi_na_index <- readRDS( file = "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/tsi_na_index.rds")

# subset tsi object with only layers containing values
lapply(seq_along(tsi_list),
       function(i) 
         tsi_list[[i]] <<- tsi_list[[i]][[ tsi_na_index[[i]] ]] 
)


# filter out reflectance values below and above 0-100 % range
# lapply(seq_along(tsi_list),
#        function(i) 
#          tsi_list[[i]] <- tsi_list[[i]] |> 
#   filter(across(1:length(names(tsi_list[[i]])), ~ . > 0),
#          across(1:length(names(tsi_list[[i]])), ~ . < 10000))
# )

3.2.1 Check NA layer occurance

# lapply(seq_along(tsi_list),
#        function(i) 
#        tsi_obs_length[[i]] - length(names(tsi_list[[i]]))
# )

lapply(seq_along(tsi_list),
       function(i) 
         (paste("In the", names(tsi_list)[[i]],"scenario", (1-round(( length(names(tsi_list[[i]]))/tsi_obs_length[[i]]), digits=4))*100, "% of layers consisted completely of NA values (n=", ( tsi_obs_length[[i]] - length(names(tsi_list[[i]]))), "out of",  tsi_obs_length[[i]], ")")) 
       
)

[[1]] [1] “In the base_case scenario 0 % of layers consisted completely of NA values (n= 0 out of 765 )”

[[2]] [1] “In the day_16_sigma_8_16_32 scenario 0.229999999999997 % of layers consisted completely of NA values (n= 2 out of 866 )”

[[3]] [1] “In the tsi_day_16_sigma_8_8_8_16_16_16_64 scenario 0 % of layers consisted completely of NA values (n= 0 out of 866 )”

[[4]] [1] “In the tsi_day_16_sigma_8_8_8_16_16_16_32_64 scenario 0 % of layers consisted completely of NA values (n= 0 out of 866 )”

#rasterVis::levelplot(tsi[[c(24)]], par.settings = viridisTheme)
#rasterVis::levelplot(tsi[[c(1,4,24,25:35)]], par.settings = viridisTheme)


for (i in seq_along(tsi_list)) {
  
  plot(tsi_stack[[i]],
       legend=T, col = viridis(n=100,option="D"),
       range=c(0,10000), 
       ext = ext(DE),
       plg=list(title= paste(names(tsi_list[i]), "\n", names(tsi_stack[[i]]))
       ))
  
   plot(vect(DE), add=TRUE)
       
}

4 generation of sample points (performed in Rstudio GUI)

4.0.1 Naming ROIs

roi_names <- c("Wesermarsch", "Greifswald", "Niederrhein",
               "Münsterland", "Westhavelland", "Eifel", "Vogelsberg",
               "Erzgebirge", "Mittelfranken", "Bayerischer_Wald", "Voralpen_Kempten", "Voralpen_Chiemsee")

4.0.2 select the ROI for each location type

# select the ROI for each location type
roi_list <- list()
lapply(1:length(roi_names), function(i)

  roi_list[[i]] <<- terra::sel(tsi_list[[2]][[c(674)]])

)
names(roi_list) <- roi_names


#save roi sample points raster stack

lapply(1:length(roi_list), function(i)

  terra::writeRaster(roi_list[[i]],
                     filename= paste0("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/",
                                      names(roi_list[i]),".tif"),
                     filetype='GTiff',
                     overwrite=TRUE)

)

4.0.3 Load the ROI for each location type

#load sampled points
roi_list <- list()
for (i in 1:length(roi_names)) {

  link <- paste0("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/",
                                      roi_names[i],".tif")
  
  roi_list[[i]] <- rast(link)
  
}
# assigen formal names to each element in list
names(roi_list) <- roi_names


# Convert each ROI to df
roi_dfs <- list()
lapply(1:length(roi_list), function(i)
  
  roi_dfs[[i]] <<- roi_list[[i]] |> as.data.frame(xy=T) |> select(x,y)  |> mutate(type=names(roi_list[i]),
                                                                                  x = as.integer(x),
                                                                                  y = as.integer(y))
)

# Combine locations into single df
 sampled_points <- data.table::rbindlist(roi_dfs)

#save sampled points
# saveRDS(sampled_points, file = "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/sampled_points_DE.rds")

# load sampled points
sampled_points <- readRDS("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/sampled_points_DE.rds")

4.1 Plots of ROIs

4.1.1 sampled location (little red highlighed areas)

knitr::include_graphics("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/Screen Shot 2022-08-28 at 11.19.41.png")

4.1.2 Detailed plots of each ROI

lapply(seq_along(roi_list),
       function(i)
         plot(tsi_stack[[i]],
              legend=T, 
              col = viridis(n=100,option="D"),
              range=c(0,10000),
              ext=ext(roi_list[[i]]),
             plg=list(title= paste(roi_names[i], "\n", names(tsi_list[[2]][[430]]))
                                 ))
)

4.1.3 compairison of TSI cases for each location

lapply(seq_along(roi_list),
       function(i)
         plot(tsi_stack,
              legend=T, 
              col = viridis(n=100,option="D"),
              ext=(ext(roi_list[[i]])),
               plg=list(title= paste(names(tsi_list[i]), "\n", roi_names[i])),
              range=c(0,10000))
       
)

For the displayed date (July 2013), all plots seem to show quite congruous results. As expected the un-interpolated time series shows more data gaps, but in general the visual correspondence of reflectance values look good.

5 TSI data manipulaiton

5.1 data manipulation and subsample creation

Due to C memory usage limit issues, the conversion from spatraster to spatial data frame had to be split into small individual processing chunks, that are then subsequently joined to form the single df. First unique start and stop limits for the processing loop are defined (as cases have heterogenious lengths depending on number of NA layers and date window size). Following this the range of the processing checks is set in a list. Then the spat raster object of each case are converted into spatial df´s, via the defined subset processing chunks. Finally, all the individual chunks of a case are joined into a single spatial df object

… works but is not a pretty workaround…

5.1.1 Create processing Chunks

tsi_df   <- list("base_case"=c(), "day_16_sigma_8_16_32"=c(), "day_16_sigma_8_8_8_16_16_16_64"=c(), "day_16_sigma_8_8_8_16_16_16_32_64"=c())

start <- list()
stop <- list()
for (i in 1:length(names(tsi_list))) {
  # create processing chunks 
  start[[i]] <- seq(from = 1, to = length(names(tsi_list[[i]]))-76, by = 76)
  stop[[i]] <- seq(from = 76, to = length(names(tsi_list[[i]])), by = 76)
  
}



# create a unique length of processing chuncks based on the number of observations in each scenario
processing_chunks <- start[[1]][c(1:length(start))] # initial list
# loop though each case
for (j in 1:length(start)) {
  # loop through each start and stop combination
  for (i in 1:length(start[[j]])) {
    
    processing_chunks[[j]][i] <- list(c(start[[j]][i], stop[[j]][i]) )
    
  }
  
}

5.1.2 Old method for spatrast to df convertion

# Convert spat raster into df, via subset processing chunks

### uncomment to execute

# OLD DF CONVERTSION METHOD
lapply(seq_along(tsi_list),
       function(i)

         for (j in c(1:length(processing_chunks[[i]]))) {

           tsi_df[[i]][[j]] <<- terra::as.data.frame(tsi_list[[i]][[c(processing_chunks[[i]][[j]][1]:
                                                                        processing_chunks[[i]][[j]][2])]],
                                                     xy = TRUE, na.rm=F) |>
             # filter out rows that only contain na
             filter_at(vars(-x,-y), any_vars(! is.na(.))) |>

             left_join(sampled_points) |>
              filter(!is.na(type))
         }
)

5.2 Spatrast to df converstion

Here the each spatrast object is subset to reduce the overall size of the created data frames. The full df are not needed, as only single locations are used for comparison between location.

# NEW DF CONVERTION METHOD
lapply(seq_along(tsi_list),
       function(i)

         for (j in c(1:#2
                     nlyr(tsi_list[[i]])
                     )) {
sub_sample <- spatSample(tsi_list[[i]][[j]], size=150000,
                  method="regular",
                  na.rm=T ,
                  as.points=T,
                  values=T,
                  as.raster=T )

#plot(sub_sample)
sub_sample_df <-terra::as.data.frame(sub_sample, geom="wkt") |>
  mutate(x = as.numeric(str_sub(geometry, 8,14)),
         y = as.numeric(str_sub(geometry, 23,29))) |>
  select(-geometry)


tsi_df[[i]][[j]] <<- sub_sample_df |> left_join(sampled_points, by = c("x", "y")) |>
              filter(!is.na(type))


 }
)

5.2.1 load output from converstion

#save tsi df
#saveRDS(tsi_df, file = "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/tsi_DE_df.rds")

# load tsi df
tsi_df <- readRDS("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/tsi_DE_df.rds")

5.2.2 Join tsi lists and reduce to single location

# create join for list into df
tsi_df_joined <- list()
 lapply(seq_along(tsi_df),
        function(i)
         tsi_df_joined[[i]] <<-  reduce(tsi_df[[i]], full_join) |>
          group_by(type) 
 )
 #save tsi df
saveRDS(tsi_df_joined, file = "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/tsi_DE_df_joined.rds")
 
 tsi_df_joined_single_location <- list()
 lapply(seq_along(tsi_df),
        function(i)
         tsi_df_joined_single_location[[i]] <<-  tsi_df_joined[[i]] |> 
          sample_n(1)
 )
 #save tsi df for single location
saveRDS(tsi_df_joined_single_location, file = "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/tsi_df_joined_single_location.rds")

5.2.3 load output from join

# load tsi df
tsi_df_joined <- readRDS("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/tsi_DE_df_joined.rds")

# load tsi df
tsi_df_joined_single_location <- readRDS("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/tsi_df_joined_single_location.rds")

5.2.4 Create long and formatted version of df for final vis

# create long version of df with unique and formatted location and date cols
tsi_df_long <- list()
lapply(seq_along(tsi_list),
       function(i) 
         tsi_df_long[[i]] <<- tsi_df_joined[[i]] |> 
         pivot_longer(cols = c(-x,-y,-type), names_to = "timestamp") |> 
         mutate(xy= paste0(x,"_",y)) |> 
         mutate(date = as.Date(timestamp, format="%Y%m%d"),
                year = lubridate::year(date),
                value = if_else(value<0,0,as.numeric(value)),
                case  = names(tsi_list[i]))
       
)


# create long version of df with unique and formatted location and date cols
tsi_df_single_location_long <- list()
lapply(seq_along(tsi_list),
       function(i) 
         tsi_df_single_location_long[[i]] <<- tsi_df_joined_single_location[[i]] |> 
         pivot_longer(cols = c(-x,-y,-type), names_to = "timestamp") |> 
         mutate(xy= paste0(x,"_",y)) |> 
         mutate(date = as.Date(timestamp, format="%Y%m%d"),
                year = lubridate::year(date),
                value = if_else(value<0,0,as.numeric(value)),
                case  = names(tsi_list[i]))
       
)

6 Visualization

6.1 Faceted time series

6.1.1 plot faceted time series 1984:2000

lapply(seq_along(tsi_list),
       function(i) 
         
         tsi_df_single_location_long[[i]] %>%
         filter(year<2001) %>%
         
         ggplot( aes(date, value, color=type, group=xy)) +
         geom_line(aes(group=type)) +
         geom_point(size=0.5) +
         theme_classic() +
         scale_colour_viridis_d(option = "D") +
         scale_x_date(date_breaks = "1 month", date_labels = "%m")+
         theme_minimal() +
         theme(axis.text.x = element_text(angle = 45,  hjust=1),
               legend.position="bottom"
         ) +
         ggtitle(cases[i]) +
         facet_wrap(~year, scales = "free_x")
       
)

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

6.1.2 plot faceted time series 2001:2021

lapply(seq_along(tsi_list),
       function(i) 
         
         tsi_df_single_location_long[[i]] %>%
         filter(year>2001) %>%
         
         ggplot( aes(date, value, color=type, group=xy)) +
         geom_line(aes(group=type)) +
         geom_point(size=0.5) +
         theme_classic() +
         scale_colour_viridis_d(option = "D") +
         theme_minimal() +
         scale_x_date(date_breaks = "1 month", date_labels = "%m")+
         theme(axis.text.x = element_text(angle = 45,  hjust=1),
               legend.position="bottom"
         ) +
         ggtitle(cases[i]) +
         facet_wrap(~year, scales = "free_x")
       
)

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

lapply(seq_along(tsi_list),
       function(i) 
         ggplot(tsi_df_single_location_long[[i]], aes(date, value)) +
         geom_smooth() +
         theme_minimal() +
         scale_x_date(date_breaks = "1 year", date_labels = "%d-%m-%y")+
         theme(axis.text.x = element_text(angle = 45 ,  hjust=1)) 
)

lapply(seq_along(tsi_list),
       function(i) 
         
         tsi_df_single_location_long[[i]] %>%
         
         
         ggplot( aes(date, value, group = year)) +
         geom_boxplot() +
         theme_minimal() +
         scale_x_date( date_labels = "%d-%m-%y")+
         theme(axis.text.x = element_text(angle = 45, hjust=1)) 
)

6.1.3 Location and Case compairison using full df

tsi_df_joined <- reduce(tsi_df_long, full_join)

## Joining, by = c("x", "y", "type", "timestamp", "value", "xy", "date", "year",
## "case")
## Joining, by = c("x", "y", "type", "timestamp", "value", "xy", "date", "year",
## "case")
## Joining, by = c("x", "y", "type", "timestamp", "value", "xy", "date", "year",
## "case")

        ggplot(tsi_df_joined, aes(x=type, y=value, group=type, fill= type, color=type)) +
  geom_boxplot(alpha=0.3) +
#  geom_jitter(alpha=0.1) +
  scale_colour_viridis_d(option = "D") +
  scale_fill_viridis_d(option = "D") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 10000)) +facet_wrap(~case)

ggplot(tsi_df_joined, aes(value, color=type, fill=type)) +
  geom_density(alpha = 0.05) +
  #geom_jitter() +
  scale_colour_viridis_d(option = "D") +
  scale_fill_viridis_d(option = "D") +
  scale_x_continuous(expand = c(0, 0)) +
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.02)) +
  theme_minimal() +
  guides(col = guide_legend(nrow = 3))+
  theme(legend.position = "bottom")  +facet_wrap(~case)

ggplot(tsi_df_joined, aes(x = value, y = type, fill = stat(x))) +
  geom_density_ridges_gradient(scale = 1.2, rel_min_height = 0.01) +
  scale_fill_viridis_c(name = "value", option = "D") +
  #scale_x_continuous(expand = c(0, 0), limits = c(0, 75)) +
  theme_minimal() +
  facet_wrap(~case)

## Picking joint bandwidth of 189
## Picking joint bandwidth of 152
## Picking joint bandwidth of 139
## Picking joint bandwidth of 139

6.2 direct location and case compairison

6.2.1 Smoothed long term trend overview

tsi_df_single_location_joined <- reduce(tsi_df_single_location_long, full_join)

## Joining, by = c("x", "y", "type", "timestamp", "value", "xy", "date", "year",
## "case")
## Joining, by = c("x", "y", "type", "timestamp", "value", "xy", "date", "year",
## "case")
## Joining, by = c("x", "y", "type", "timestamp", "value", "xy", "date", "year",
## "case")

ggplot(tsi_df_single_location_joined, aes(date, value, color=case, group = xy)) +
  geom_line(aes(group=case), alpha=0.5) +
  geom_point(size=0.5, alpha=0.9) +
  geom_smooth(color="black") +
  theme_minimal() +
  scale_colour_viridis_d(option = "D") +
  scale_x_date(date_breaks = "4 year", date_labels = "%Y")+
  theme(axis.text.x = element_text(angle = 45,  hjust=1),
        legend.position="bottom") +
  #  ggtitle(cases[i]) +
  facet_grid(rows = vars(type), cols = vars(case), scales="free_y") 

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

#facet_grid( ~year)
#facet_grid(vars(drv), vars(cyl))

6.2.2 compairson of Base case values and interpolated TSI outputs

year_chunks <- round(seq(1984, 2021, length.out = 5))
lapply(seq_along(c(1:4)),
       function(i) 

tsi_df_single_location_joined %>%
  filter(year > year_chunks[i],
         year < year_chunks[i+1],
         case != "base_case") %>%
  
  ggplot(aes(date, value, color=case, group = xy)) +
  geom_line(aes(group=case)) +
  geom_point(size=0.5) +
  geom_point(data=tsi_df_single_location_joined |> filter(case =="base_case",
                                                 year > year_chunks[i],
                                                 year < year_chunks[i+1]),
             aes(date, value), 
             color ="red", alpha=0.5) +
  theme_minimal() +
  scale_colour_viridis_d(option = "D") +
  scale_x_date(date_breaks = "1 month", date_labels = "%m")+
  theme(axis.text.x = element_text(angle = 45,  hjust=1),
        legend.position="bottom") +
  facet_grid(rows = vars(type), cols = vars(year), scales="free_x") 

#  facet_wrap(~year, scales = "free_x") 

)

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

6.2.3 plot time series with formatted dates

ggplot(tsi_df_long[[i]], aes(date, value, group=xy, color=xy)) +
  geom_line() +
  theme_classic() +
  scale_colour_viridis_d(option = "D") +
  scale_x_date(date_breaks = "1 month", date_labels = "%m-%y")+
  theme(axis.text.x = element_text(angle = 45,  hjust=1),
        legend.position="none")

ggplot(tsi_df_long[[i]], aes(date, value)) +
  geom_smooth() +
  theme_minimal() +
  scale_x_date(date_breaks = "1 year", date_labels = "%d-%m-%y")+
  theme(axis.text.x = element_text(angle = 45 ,  hjust=1)) 

ggplot(tsi_df_long[[i]], aes(date, value, group = timestamp)) +
  geom_boxplot() +
  theme_minimal() +
  scale_x_date( date_labels = "%d-%m-%y")+
  theme(axis.text.x = element_text(angle = 45, hjust=1))